##########################################################################################

library(data.table)
library(optparse)
library(ggplot2)
library(Seurat)
library(ggsci)
library(ArchR)
library(ggplotify)

##########################################################################################
option_list <- list(
    make_option(c("--rna_file"), type = "character"),
    make_option(c("--atac_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 单细胞表达文件
    rna_file <- "~/20231121_singleMuti/results/qc_atac/testis_combined.annotationCellType.qc.Rdata"

    ## atac文件
    atac_file <- "~/20231121_singleMuti/results/qc_atac/testis_combined_peak.combineRNA.qc.Rdata"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/celltype_plot/tmp"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

rna_file <- opt$rna_file
atac_file <- opt$atac_file
out_path <- opt$out_path

dir.create( paste0(out_path , "/markerGenes") , recursive = T)

###########################################################################################

a <- load(rna_file)
DefaultAssay(scrnat) <- "MAGIC_RNA"
## scrnat

b <- load(atac_file)
## testis_combined_peak_combineRNA

###########################################################################################

marker_genes <- c("GFRA1", "UTF1", "KIT", "DMRT1", "STRA8",
"TNP1", "PRM2", "NOTCH3", "OVOL1", "OVOL2",
"DLK1", "CD68", "CD14", "MYH11", "ACTA2",
"SPO11", "SYCP3", "AMH", "SOX9", "VWF",
"PECAM1", "NME8", "C9orf116", "ACR", "TXNDC8",
"TXNDC2", "NKG7", "FGFBP2","SYCP3")

## 细胞颜色
use_colors <- c(pal_npg("nrc")(10) , pal_jco("default")(6))
#names(use_colors) <- unique(scrnat$cell_type)
names(use_colors) <-c( "Myoid cells","Leydig cells" ,          
"Endothelial cells","Zygotene",         
"Round&ElongateS.tids","Patchytene",
"SSC","Sperm" ,                 
"Diplotene","Early stage of spermatids", 
"Leptotene","Sertoli cells",            
"Macrophages","Differenting&Differented SPG",
"Pericytes","NKT cells")

## 样本类型颜色
col_sample <- c(
    rgb(red=212,green=31,blue=37,alpha=255,max=255) ,
    rgb(red=31,green=137,blue=66,alpha=255,max=255) ,
    rgb(red=39,green=45,blue=106,alpha=255,max=255) 
    )
names(col_sample) <- c("testis01" , "testis02" , "testis03")

###########################################################################################
## 画聚类图 ##
## 细胞类型
#scrnat <- RunUMAP(scrnat, dims = 1:10)

out_file <- paste0( out_path , "/magic3_cluster.qc.cell_type.pdf" ) 
p <- DimPlot(scrnat, reduction = "umap", label = TRUE, group.by = 'cell_type' , cols = use_colors ) + theme(legend.position = 'none') 
ggsave( out_file , p ,width = 6 , height = 6)
out_file <- paste0( out_path , "/magic3_cluster.qc.cell_type.legend.pdf" ) 
p <- DimPlot(scrnat, reduction = "umap", label = TRUE, group.by = 'cell_type' , cols = use_colors  ) 
ggsave( out_file , p ,width = 6 , height = 6)

## 去除umap不合适的图
umap_type = scrnat@reductions$umap@cell.embeddings %>% as.data.frame() %>% cbind(type = scrnat@meta.data$cell_type)

## 不同来源
out_file <- paste0( out_path , "/magic3_cluster.qc.sample.pdf" ) 
p <- DimPlot(scrnat, reduction = "umap", label = FALSE, group.by = 'batch' , cols = col_sample  ) + theme(legend.position = 'none') 
ggsave( out_file , p ,width = 6 , height = 6)
out_file <- paste0( out_path , "/magic3_cluster.qc.sample.legend.pdf" ) 
p <- DimPlot(scrnat, reduction = "umap", label = FALSE, group.by = 'batch' , cols = col_sample  ) 
ggsave( out_file , p ,width = 6 , height = 6)


## marker基因 ##
for(gene in marker_genes){
    out_file <- paste0( out_path , "/markerGenes/" , gene , "_cluster_marker.qc.pdf" ) 
    p <- FeaturePlot( scrnat , features = gene, pt.size = 0.2 , cols= c("grey","red") ) +
        theme(
            legend.position = 'none',
            legend.title = element_blank() ,
            panel.grid.major=element_blank(),
            panel.grid.minor=element_blank(),
            panel.background = element_blank(),
            panel.border = element_rect(color = "black", size = 0.5),
            plot.title = element_text(size = 12,color="black",face='bold'),
            axis.line = element_blank(),  # 隐藏 x 轴横线
            axis.text = element_blank(),  # 隐藏 x 轴文本
            axis.title = element_blank(),  # 隐藏 x 轴标题
            axis.ticks.length = unit(0, "cm") ## x轴坐标不显示
        ) 
    ggsave( out_file , p ,width = 3 , height = 3)
}

out_file <- paste0( out_path , "/markerGenes/legend.pdf" ) 
p <- FeaturePlot( scrnat , features = gene, pt.size = 0.2 , cols= c("grey","red") ) +
    theme(
        legend.position = 'right',
        legend.title = element_blank() ,
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(color = "black", size = 0.5),
        plot.title = element_text(size = 12,color="black",face='bold'),
        axis.line = element_blank(),  # 隐藏 x 轴横线
        axis.text = element_blank(),  # 隐藏 x 轴文本
        axis.title = element_blank(),  # 隐藏 x 轴标题
        axis.ticks.length = unit(0, "cm") ## x轴坐标不显示
    ) 
ggsave( out_file , p ,width = 3 , height = 3)

#精原细胞(UTF1) —Cell Research. 2018
#正在分化的精原细胞(KIT)—Cell stem cell. 2018
#分化完成的精原细胞(STRA8)—Cell stem cell. 2018
#细线期精母细胞(SPO11)—Cell stem cell. 2018
#偶线期精母细胞(SYCP3)—Cell Reports. 2018
#粗线期精母细胞(OVOL2)—Cell stem cell. 2018
#双线期精母细胞(NME8)—Cell stem cell. 2018
#圆形精子和长形精子细胞(TXNDC8)		—Cell stem cell. 2018
#精子(TNP1、PRM1)—Cell stem cell. 2018
#Sertoli细胞(AMH)—Cell stem cell. 2018
#Leydig细胞(DLK1)—Cell Research. 2018
#肌样细胞(MYH11)—Cell stem cell. 2018
#周细胞(NOTCH3)—Human Molecular Genetics.2022
#巨噬细胞(CD14)—Cell Research. 2018
#内皮细胞(VWF)—Cell Research. 2018
#NKT细胞(NKG7、FGFBP2)

###########################################################################################
## ATAC
projHeme5 <- testis_combined_peak_combineRNA

## atac的聚类图,标记细胞类型
p <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "cell_type", embedding = "UMAP", plotAs="points" , pal=use_colors)
#p1 <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "Clusters", embedding = "UMAP")
p2 <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "Sample", embedding = "UMAP", plotAs="points" , pal=col_sample)

p <- as.ggplot(p)
out_file <- paste0( out_path , "/plotEmbedding.atac.qc.cell_type.pdf")
ggsave( out_file , p ,width = 7 , height = 7)
p2 <- as.ggplot(p2)
out_file <- paste0( out_path , "/plotEmbedding.atac.qc.sample.pdf")
ggsave( out_file , p2 ,width = 7 , height = 7)
